{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1.7 Complex-valued waves\n",
    "====\n",
    "\n",
    "\n",
    "In NGSolve finite element spaces can be built and linear systems can be solved over the complex field. This tutorial shows how to compute the solution of the Helmholtz equation with impedance boundary conditions in complex arithmetic. The boundary value problem is to find $u$ satisfying \n",
    "\n",
    "$$\n",
    "-\\Delta u - \\omega^2 u = f\\qquad \\text{ in } \\Omega\n",
    "$$\n",
    "\n",
    "together with the impedance (outgoing) boundary condition\n",
    "\n",
    "$$\n",
    "\\frac{\\partial u }{ \\partial n} - i \\omega u = 0 \n",
    "\\quad \\text{ on } \\partial \\Omega\n",
    "$$\n",
    "\n",
    "where $i =$ `1j` is the imaginary unit.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.geom2d import SplineGeometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Geometry \n",
    "geo = SplineGeometry()\n",
    "geo.AddCircle((0.5, 0.5), 0.8,  bc=\"outer\")\n",
    "geo.AddRectangle((0.7, 0.3), (0.75, 0.7),\n",
    "                 leftdomain=0, rightdomain=1, bc=\"scat\")\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.05))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Declare  a complex finite element space "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=5, complex=True)\n",
    "u, v = fes.TnT()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Wavenumber & source\n",
    "omega = 100\n",
    "pulse = 1e3*exp(-(100**2)*((x-0.5)*(x-0.5) + (y-0.5)*(y-0.5)))\n",
    "Draw(pulse, mesh, 'pulse')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Forming the system\n",
    "\n",
    "The weak form for $u \\in H^1$:\n",
    "\n",
    "$$\n",
    "\\int_\\Omega\\big[ \\nabla u \\cdot \\nabla \\bar v - \\omega^2 u \\bar v \\big]\n",
    "\\, dx - \n",
    "i \\,\\omega\\, \\int_{\\partial \\Omega} u \\bar v \\, ds = \\int_{\\Omega} f \\bar v\n",
    "$$\n",
    "\n",
    "for all $v$ in $H^1$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Forms\n",
    "a = BilinearForm(fes)\n",
    "a += grad(u)*grad(v)*dx - omega**2*u*v*dx\n",
    "a += -omega*1j*u*v * ds(\"outer\")\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(fes)\n",
    "f += pulse * v * dx\n",
    "f.Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes, name=\"u\")\n",
    "gfu.vec.data = a.mat.Inverse() * f.vec\n",
    "Draw(gfu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Explore the GUI's menu options in `Visual` tab:\n",
    "\n",
    "- Increase subdivions\n",
    "- Real and imaginary parts\n",
    "- View absolute value\n",
    "- Turn off Autoscale\n",
    "- Turn on Deformation\n",
    "- Turn on Periodic Animation\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
